% Define the stability functions for s=1,2,3,4
R1 = @(z) 1 + z;
R2 = @(z) 1 + z + z.^2/2;
R3 = @(z) 1 + z + z.^2/2 + z.^3/6;
R4 = @(z) 1 + z + z.^2/2 + z.^3/6 + z.^4/24;
R5 = @(z) 1 + z + z.^2/2 + z.^3/6 + z.^4/24 + z.^5/120;
R6 = @(z) 1 + z + z.^2/2 + z.^3/6 + z.^4/24 + z.^5/120 + z.^6/720;
% Create a grid in the complex plane
x = linspace(-4, 2, 500);
y = linspace(-4, 4, 500);
[X, Y] = meshgrid(x, y);
Z = X + 1i*Y;

% Compute |R(z)| for each method
absR1 = abs(R1(Z));
absR2 = abs(R2(Z));
absR3 = abs(R3(Z));
absR4 = abs(R4(Z));
absR5 = abs(R5(Z));
absR6 = abs(R6(Z));

% Plot the stability regions (|R(z)| = 1)
figure;
hold on;
contour(X, Y, absR1, [1, 1], 'r', 'LineWidth', 1, 'LineStyle', ':');   % s=1
contour(X, Y, absR2, [1, 1], 'g', 'LineWidth', 1, 'LineStyle', '-.');   % s=2
contour(X, Y, absR3, [1, 1], 'b', 'LineWidth', 1);   % s=3
contour(X, Y, absR4, [1, 1], 'm', 'LineWidth', 1, 'LineStyle', '-');   % s=4
contour(X, Y, absR5, [1, 1], 'y', 'LineWidth', 1, 'LineStyle', '--');   % s=5
contour(X, Y, absR6, [1, 1], 'c', 'LineWidth', 1);   % s=6
hold off;

% Add labels and legend
grid on;
xlabel('Re(z)');
ylabel('Im(z)');
title('Stability Regions of Explicit RK Methods (s=1,2,3,4,5,6)');
legend('s=1', 's=2', 's=3', 's=4','s=5','s=6', 'Location', 'southwest');
axis equal;
axis([-4 2 -4 4]);